Prepare the R environment

We load the packages needed for this practice.

## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo

The first step is to load the data from the the CABIN Canadian Aquatic Biomonitoring Network in the Hamilton Harbour. We then want to project in the NAD83 / Ontario MNR Lambert (more accurate and suitable for the study area). If you want to read more about it: http://epsg.io/3161

Load data

The .rds can be load in the R environment using the function readRDS(). The file contains an R object: the CABIN data for the Hamilton Harbour.

Cloud variogram

Directional variogram

To convert the variogram into a directional one, we have to set the argument alpha. We pass a vector of angles through the alpha parameters: c(0, 45, 90, 135).

Empirical variogram

Emprical sometimes referred to sample variogram.

Step 2 - Cloud variogram

Construct a cloud variogram with a maximum distance of 1500m.

Step 3 - Check for outliers

Select point pairs by digitizing a region with the mouse (left mouse button adds a point, right mouse button ends)

Ids of the row for the pair of points

##   head tail
## 1    4    6
## 2    5    6

Step 4 - Remove outlier

Residual empirical variogram

It is also possible to study the autocorrelation structure in the data after accounting for one or more covariates.

In the example below, we evaluate the autocorrelation structure after accounting for sampling effort.

## Warning in variogram.default(y, locations, X, trend.beta = beta, grid = grid, :
## linear model has singular covariance matrix

Exercice

  1. Build a model variogram using the empirical variogram
  • This means using another model variogram fitting method
  1. Build a model variogram with another variable
  • This may require using another model variogram
  1. Build a model variogram on your own data
  • Think of the question you want to answer
  • The size of your data will influence the choice of model fitting procedure

Kriging

Using the model variogram defined above we can now build kriging maps

Step 1 - Define a grid

## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
## class      : RasterLayer 
## dimensions : 82, 148, 12136  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 1341058, 1355858, 11862377, 11870577  (xmin, xmax, ymin, ymax)
## crs        : +proj=lcc +lat_1=44.5 +lat_2=53.5 +lat_0=0 +lon_0=-85 +x_0=930000 +y_0=6430000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

As for inverse distance weighting, we build a grid of \(100 \times 100\) m.

Ordinary kriging

Draw the map

Exercice

  1. Build a kriging maps with another variable
  • This means using another model variogram
  1. Perform universal kriging
  2. Build a kriging maps on your own data
  • Think of the question you want to answer
  • The size of your data will influence the choice of model fitting procedure